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Abstract 

The space-time dynamics of rigid inhomogeneities (inclusions) free to move in a ran- 
domly fluctuating fluid bio-membrane is derived and numerically simulated as a function 
of the membrane shape changes. Both vertically placed (embedded) inclusions and hori- 
zontally placed (surface) inclusions are considered. The energetics of the membrane, as a 
two-dimensional (2D) meso-scale continuum sheet, is described by the Canham-Helfrich 
Hamiltonian, with the membrane height function treated as a stochastic process. The 
diffusion parameter of this process acts as the link coupling the membrane shape fluc- 
tuations to the kinematics of the inclusions. The latter is described via Ito stochastic 
differential equation. In addition to stochastic forces, the inclusions also experience 
membrane-induced deterministic forces. Our aim is to simulate the diffusion-driven ag- 
gregation of inclusions and show how the external inclusions arrive at the sites of the 
embedded inclusions. The model has potential use in such emerging fields as designing 
a targeted drug delivery system. 

PACS 87.20- Membrane biophysics. 

PACS 34.20- Interatomic and intermolecular potential. 

PACS 87.22BT- Membrane and subcellular physics. 
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Amphiphilic molecules, such as lipids and proteins, can self-assemble themselves into a variety 
of exotic structures in an aqueous environment |T|]. Computer-based simulation of the dynamics 
of these structures forms an interesting research area in the computational statistical mechanics 
of bio- material systems. One such structure is the phospholipid bilayer which represents the 
generic structure of all bio-membrane systems, both natural and artificial. These membranes 
can have thicknesses of only few nano-metres but linear sizes of up to tens of micro-metres, 
and can therefore be regarded as highly flexible, fluid-like, 2D continuum sheets embedded in 
a three-dimensional space. 

Thermal fluctuations can induce shape fluctuations and shape transformations in mem- 
branes. For example, in the so-called budding transition a spherical vesicle transforms into 
a prolate ellipsoid as the temperature increases. There is also the possibility that the spherical 
geometry becomes oblate, producing a shape similar to the biconcave rest shape of a red blood 
cell. 

Bio-membranes regulate the recognition and transport processes to and from the interior 
of the cells, as well as between the cells, forming a barrier which all external particles arriving 
at the cell must cross. They contain a variety of integral (embedded) inhomogeneities, such as 
proteins and other macromolecules , that penetrate the thickness of the membrane and act 
as transport channels. We shall refer to these as the M-type inclusions. These inclusions are 
mobile and can freely diffuse across the membrane. Their presence, however, force the bilayer 
to adjust its thickness locally so as to match the thickness of the hydrophobic region of these 
inclusions ^, causing local deformations in the membrane geometry. The perturbations 
produced in the membrane shape due to these local deformations give rise to both short- 
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range and long-range membrane-induced indirect forces between the inclusions. These forces 
act in addition to direct Van der Waals and electrostatic forces between the inclusions. The 
Long-range forces originate from the perturbations associated with the long wavelength shape 
fluctuations j^], whereas the short-range forces are associated with the local deformations in the 
immediate vicinity of the inclusions . These membrane- induced forces between the inclusions 
play a far more significant role, viz the direct molecular interactions, when the length scales 
involved are comparable to the size of the membrane. In addition to this mode of deformation, 
a membrane can also deform due to the tension at the amphiphilic molecules-water interfaces. 
This tension results in a change in the overall surface area of the membrane. A third mode 
of deformation also exists and this is associated with the bending elastic-curvature property 
of the membrane which distinguishes it from a sheet of simple fluid dominated by surface 
tension f^. Accordingly, two models to study the inclusion-induced local deformations have 
been developed. In the first model, the membrane energy is taken to consist of a contribution 
from the molecular expansion/compression due to the change in the thickness at the inclusion 
boundary, and also a contribution from the overall change in the surface area. Using this 
model, it is shown that |]^, |^, |lT| the inclusion- induced deformations cause exponential 



decays in the thickness of the membrane, extending from the inclusion-imposed value to the 
equilibrium thickness value, as shown schematically in Fig.l for two rod-like inclusions. In the 



second model [^, |I2[, the contribution of the membrane bending property is taken into account 
in the energy term, and it is found that this significantly affects the perturbation profile at 
the inclusion boundary as well as modifying the membrane- induced interactions. Evidently, 
an object supported by surface tension would have a different dynamics than one supported 
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by the bending elasticity of the surface [|T^. 

In addition to the M-type inclusions, membranes can also carry inclusions that lie on their 



surfaces as shown schematically in Fig. 2. These surface inclusions can represent objects 
that have arrived at the membrane from the outside and are therefore referred to as external 
inclusions. We shall refer to these as the S-type inclusions. 

At the meso-scale, i.e. when the detailed molecular architecture of the membrane can be 
subsumed into a background 2D sheet, the free elastic energy of a symmetric membrane is 



described by the Canham-Helfrich Hamiltonian |T^, [IB 



n = J <fa^{ao + 2kH'^ + RK} , (1) 



where 



K = detiK.j) = —- 
Kin 

1„ . 1/1 1 
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i/ = -Ti.(A-,) = -^^ + ^). (2) 



are respectively the Gaussian and mean curvatures of the sheet, Ri and R2 are the two principle 
radii of curvature of the sheet, cxo is the surface tension, n is the bending rigidity, R is the 
Gaussian rigidity, g is determinant of the metric tensor and a = (cri,o"2) is the 2D local 
coordinate on the sheet as opposed to the coordinates on the embedding space. The last term 
in (|1|) is, by Gauss-Bonett theorem, an invariant for closed surfaces implying that the dynamics 
of a membrane is not influenced by this term if its topology remains fixed. In what follows, 
we concentrate on membranes with fixed topology and drop this term. We then have 

n = J rfV^{(7o + 2kH^}. (3) 



The study of a membrane whose free energy is described by is facihtated by considering it 
to be nearly flat, i.e. its thickness to be much smaller than its linear size L. This is indeed 
what we mean by a meso-scale model of a membrane. We therefore take the membrane to be 
almost parallel to the (xi,a;2) plane, regarded as the reference plane. The position of a point 
on the membrane can then be described by a single- valued function h{xi^X2) representing the 
height of that point. This simplification is achieved by writing the Hamiltonian (^) in the 
Monge representation which gives for the mean curvature 



2H = -g-'^'^ldlh^l + {d2hf) + dlh{l + {d^hf) - 2dihd2hd^d2h] , (4) 

where di = i = 1,2. We assume that the area of the membrane can fluctuate without 
constraint by setting ao = in (|^). Consequently, using (p, the Hamiltonian to leading 
order in derivatives of h becomes 

no = ^J dM^'hf. (5) 

This is the Canham-Helfrich Hamiltonian in Monge representation, expressed in terms of the 
height function of the membrane. It is the expression that we employ to describe the energetics 
of the membrane. 

Employing a statistical mechanics based on (|^) only, and ignoring the contributions from 
the expansion/compression and interfacial energies, the potential energy function 

\2A^ 

^MMiRij) = -^bT-—^ , (6) 

was constructed to describe the membrane-induced temperature- dependent long-range 
forces between a pair of disk shape M-type inclusions that can freely tilt with respect to each 



other. Another function was also constructed for long-range interaction between two S-type 
inclusions [0] 

VIsiR^v 0., 9,) = -ksT^^ cos-'m + 0,)] , (7) 

where A = tttq is the area of an M-type inclusion of radius ro, ks is the Boltzmann constant, 
Rij is distance between the centres of mass of two inclusions i and j, and Lj are the lengths 
of two S-type inclusions making the angles 9i and 9j respectively with the line joining their 
centres of mass (see Fig. 2) and T is the membrane temperature. It is evident that both of 
these membrane-induced potentials are attractive and fall off as with the distance. These 
expressions are derived for rod-like inclusions that are assumed to be much more rigid than 
the ambient membrane so that these inclusions can not move coherently with the membrane. 
The only degrees of freedom for the rods are rigid translations and rotations while they remain 
attached to the membrane. 

So far, the modelling of bio-membrane dynamics decorated with inclusions has been mainly 
concerned with constructing potential energy functions such as those given in (|^) and (0). 
An interesting problem, however, would be to use this information to simulate the space- 
time behaviour of inclusions in a membrane described by (^) and undergoing stochastic shape 
fluctuations. This is the problem that we address in this paper. This type of simulation can 
establish a direct link between the randomly changing membrane shape on the one hand and the 
inclusion dynamics on the other. In such a simulation, the thermodynamic phase behaviour 
of inclusions, such as their temperature-dependent aggregation, can be directly computed. 
This phase behaviour plays a crucial role in the functional specialisation of a membrane 0. 
Furthermore, information on the capture rate of the S-type inclusions, which could represent 
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external drug particles, at the sites of the M-type inclusions can be obtained as a function of the 
changes in the environmental variables such as the ambient temperature. This type of meso- 
scale simulation when coupled with the Molecular Dynamics (MD) simulation of membrane 



patches near the inclusions at the nano-scale |T^, can produce a seamless multi-scale model of 



the entire environment for many bio-molecular processes, starting with the arrival of external 
inclusions at the cell, their diffusion in the membrane, and finally their molecular docking at 
the site of the embedded inclusions. 

To proceed, let us consider a 2D bio-membrane described by (|^) containing both the M- 
type and the S-type inclusions. To make the membrane a stochastically fluctuating medium, 
we treat the height function in as a stochastic Wiener process with a Gaussian distribution, 
whose mean and variance can be written as 

(/i(xi,X2;t)) = 0, (8) 

{h{xi,X2;t) h{xi,X2;t)) = {[h{xi,X2;t) - {h{xi,X2;t))]'^) = 2Dt (9) 

where D is the diffusion constant associated with the height fluctuations at the local position 
{xi,X2) and represents the measure with which random fluctuations propagates in the local 
geometry. Such random height changes would cause a roughening of the membrane surface on 
molecular scales, and this has been observed in NMR experiments . 



Assuming that this is the only stochastic process present in the membrane, it is then 
reasonable to suppose that this stochastic dynamics is communicated to the inclusions as well, 
and that their ensuing random motions are contingent only on these fluctuations. This implies 
that the mathematical point representing the centre of mass of an inclusion coinciding with 
the membrane point {xi,X2) would also experience the same fluctuations and would diffuse 
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with the same diffusion constant. To derive an expression for D, based on we start with 
the static height-height correlation function obtained from (j^). This is given by [Q, ^ 

(Mq; 0)/i*(q'; 0)) = ^(27r)^5(q - q') , (10) 

where (■ ■ ■) is the thermal averaging with respect to the Boltzmann weight, exp {— Tio / ksT) , 
and q is the wave vector of magnitude q. The corresponding dynamic correlation function can 
be obtained by writing 

/i(q; t) = /i(q; 0)e~^"('i)* , (11) 

giving 

{h{q;t)h*{q!;t)) = ^e-'^^^^^\2nf5ici-c^) , (12) 

where the damping rate, 'Jo{q), reflecting the long-range character of the hydrodynamic damp- 
ing, is defined as 

1o{q) = Kq^/^V, (13) 

and 7] denotes the coefficient of viscosity of the fluid membrane. In real space, (p!2[) transforms 
00 to 



{h{xi, X2; t)h{xi,X2; t)) = -^L^e'^^"* , (14) 

AlTK 

where L is the length of the membrane. This is the equal-time correlation function for mem- 
brane fluctuations. A similar model of an active fluctuating membrane in which the vertical 
displacements of the membrane satisfy a Langevin equation in the q space has also been pro- 
posed and is its shown that a term similar to the static version of (|I2p contributes to 



the correlation function which also contains a contribution from non-equilibrium fluctuations. 
The latter is in the form of a q~^ term which dominates at long distances. Comparison of (^ 



and ( |T^ ) yields the desired result 



2g-27ot 



D = ^^^^^ . (15) 

It should be emphasised that the association of a diffusive process with the membrane height 
function, and the resulting diffusion constant, is not analogous to the usual model of a diffusion 
process in which, for example, a particle diffuses through a medium, such as a fluid. Rather, 
what is suggested here is that the magnitude of a mathematical function representing the 
height of a mathematical point in the membrane is subject to random stochastic variations, 
and the diffusion constant is a measure of this variation. 

When the M-type inclusions are present they produce exponentially decaying local defor- 
mations in the membrane geometry (see Fig.l). Correspondingly, the correlation function can 
be modified by a multiplicative exponential factor to 

{h{xi, X2\ t)h{xi, X2\ t))Ni = e~''°/^{h{xi, Xa; t)h{xi, Xa; t)) , (16) 

where tq is the radius of an M-type inclusion, R + vq is the radius of the circular region 
around the inclusion with its centre coinciding with that of the inclusion, and NI stands for 
near inclusion. It is evident that outside this region the exponential decay of the profile is 
negligible. Accordingly, within this circular region of radius R + tq, the diffusion constant is 
also modified to 

= De-^»/^. (17) 

This equation implies that when the centre of mass an S-type inclusion enters a circular region 
of radius R + tq its diffusion coefficient goes over to Dm and progressively approaches zero as 
the boundary of an M-type inclusion is approached. We can ascertain, as a first approximation, 
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that this is how an M-type inclusion interacts with an S-type inclusion. 

In our simulations, the equations motion of both the S-type and the M-type inclusions are 



represented by the differential equation of the the Ito stochastic calculus 

dr{t) = A[r(t), t] dt + D^''^d}N{t). (18) 

This equation describes the stochastic trajectory, r(t), of the centres of mass of the inclusions 
in terms of a dynamical variable of the inclusions, A[r(t),t], which is referred to as the drift 
velocity, and a term, dW{t), which is a given Gaussian Wiener process with the mean and 
variance given by 

(rfW(t)) = 

(19) 

{dWi{t)dWj{t)) = 26ijdt. 

Equation (piSf ) applies to each dimension of the motion. The Ito equation predicts the increment 
in position, i.e. dr{t) = r{t + dt) — r(t), for a meso-scale time interval dt as a combination of 
a deterministic drift part, represented by A[r{t),t], and a stochastic diffusive part represented 
by D^/'^dW{t) and superimposed on this drift part. This equation resembles the 'position' 
Langevin equation describing the Brownian motion of a particle [^. The position Langevin 



equation corresponds to the long-time (diffusive time) configurational dynamics of a stochastic 
particle in which its momentum coordinates are in thermal equilibrium and hence have been 
removed from the equations of motion. Since we are interested in diffusive time scales as well, 
we can re-write (|18|) as 

dr(t) = 7^F(t) dt + D^/^dW(t) , (20) 
11 



where F(t) is the instantaneous systematic force experienced by the i-th inclusion and is 
obtained from the inter-inclusion potentials, given in (g) and (J^), according to 

F. = -EVR,m.)- (21) 

j>i 

We implemented ( pO|) for our 2D simulations according to the iterative scheme 



X{t + dt) = X{t) + -^Fx{t) dt + V2Ddt R% 



(22) 
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Y{t + dt) = Y{t) + Fy(t) dt + V2Ddt B2 

kbT 

where Rx and Ry are standard random Gaussian variables chosen separately and indepen- 
dently for each inclusion according to the procedure given in and Fx,Fy are the X 
and Y components of the force F. For the S-type inclusions, we treated the angles in (0) as 
independent stochastic variables described by 

D 



e(t + dt) = eit) + -Tit) dt + -V2Ddt , (23) 

KbTL^ L 

where r is the torque experienced by an S-type inclusion and is given by 

_ dV'^{Rij,9i,9j) 

h 96, ' ^^^^ 

j>i * 

and 9'^ is the angular counterpart of Rx and Ry- 

In the numerical simulations, recently reported in their broad outlines ||25[ , we use a square 
membrane with L = AOfim on its side. The other parameters used were set at k = 10~^^ J 
and 7] = 10~^ J sec m^'^ [Q. These values correspond to the condition in which the bending 
mode of the membrane is important. From these data the damping coefficient, 70, in the real 
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space of the membrane, can be obtained from (pTSl). The simulation temperature was set at 
T = 300° K, and the correlation (delay) time, t, over which the diffusion coefficient in (|T3|) 
was calculated, was set at t = 10~^sec. These data gave D = 2.6 x 10^^ m^sec"^, which is in 
close agreement with the value of ~ 4.4 x lO^^m^sec"^ obtained at the molecular level via 
an MD simulation of a fully hydrated phospholipid dipalmitoylphosphatidylcholine (DPPC) 
bilayer diffusing in the ^-direction pGf. To justify our choice of the correlation time, t = 10~^ 
sec, we recall that the time scale of a stochastic particle, t^, of mass m is usually determined 
from the relation 

rriD , , 

to = ^ (25) 

Since to is normally of the order of 10~^sec, then for the criterion of long-time dynamics, 
employed in our model (cf (pO])), to be justified the correlation (diffusive) time scale, t, in (|T5|) 



has to satisfy the condition 



t>ti). (26) 

For our calculated value of D and our choice of the inclusion mass m = lO^^^/xg, corresponding 
to an inclusion of length Lj = 0.1/im, we obtained a value of to = 0.6 x lO^^sec, showing that 



our choice of the correlation time was appropriate to satisfy the condition in ([26| ) . The radius 
of an M-type inclusion was set at tq = O.Olyum, and the inclusions were all equal in length. 

The stochastic trajectories of the inclusions were obtained in a set of five simulations. The 
simulation time step, dt, in (|2^) was set at dt = lO^^sec, and each simulation was performed 
for 4 X 10^ time steps, i.e for a mesoscopic interval of 4000yusec. The total number of inclusions 
considered was 36, consisting of 13 S-type and 23 M-type. 
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In the first simulation, we computed the random motions of the S-type inclusions in a 
membrane devoid of the M-type inclusions. This was done in order to observe the details 
of the drift-diffusion motion over mesoscopic time scales. Figure 3 shows the stochastic X-Y 
trajectories of a sample of 4 S-type inclusions plotted on a micron scale up to the end of the 
simulation time. In addition to the drift motions, represented by the second terms in (p2|), 
the random Brownian-type variations, emanating from the membrane shape fluctuations, are 
superimposed on this drift motion and are clearly visible over the mesoscopic length and time 
scales. Figure 4 shows the snapshots of a small patch of the membrane with both the S-type 
(white spheres) and the M-type (black spheres) inclusions. In this, and subsequent figures, the 
solid spheres refer to the centres of mass of the rod-like inclusions. In the initial state the outer 
M-type inclusions were regularly positioned, whereas the inner ones were randomly distributed. 
The S-type inclusions were all distributed completely at random. Figures 4a to 4c refer to the 
simulation in which the M-type inclusions were pinned to the membrane, i.e. were static, and 
only the S-type inclusions were mobile, and figures 4d to 4f refer to the simulation in which 
both the M-type and the S-type inclusions were mobile. The initial states in both simulations, 
figures 4a and 4d, were the same. The snapshots were obtained from dynamic simulations, 
akin to a MD simulation, covering the entire simulation time interval. These snapshots were 
recorded at 2 x 10~^ sec intervals, with figures 4c and 4f referring to the final states reached 
at the conclusion of the simulations after 4 x 10^ time steps. The animation of a complete 
run showed clearly the stochastic motions of the inclusions, and how the S-type inclusions 
approached the M-type inclusions and were captured at the site of the M-type inclusions. An 
examination of figure 4 shows that for the case of dynamic M-type inclusions, a larger number 
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of the S-type inclusions were captured at the M-type inclusion sites, i.e. the number was some 
4 times higher than in the static case at the same temperature. We adopted the method of 
counting an S-type inclusion as a captured inclusion when its centre of mass coincided with 
that of an M-type inclusion. The numerical algorithm then transformed the colour code of that 
S-type inclusion from white to black. Figures 4e and 4d also show some diffusion-driven local 
aggregation of the M-type inclusions. The capture of the S-type inclusions can be viewed as 
the first stage in the molecular docking process which will eventually transfer these inclusions 
into the interior of the cell. 

To examine the membrane response to temperature changes, two of the simulations were 
performed at different temperatures. Figure 5 shows the snapshots of these simulations at 
T = 100°K (a to c) and at T = 350°K (d to f). Figures 5d to 5f clearly show that both the 
number of captured inclusions and the aggregation of the M-type inclusions were affected by 
these temperature differences, as can be seen by comparing figures 5c and 5f. 

To sum up, although many dynamical aspects of membrane-like surfaces have been ad- 



dressed in the past [28|, it is only relatively recently that attention has focused on the dy- 
namics of membranes with inclusions. To our knowledge no computer-based simulation of this 
dynamics has been reported so far. In this paper we constructed a meso-scale model of a 
generic bio-membrane based on the Canham-Helfrich curvature-energy formalism. We treated 
the height function of the membrane as a stochastic Wiener process whose correlation function 
provides the relevant diffusion constant describing the membrane fluctuations. Two types of 
inclusions, one mimicking the internal embedded type and the other the external floating type, 
are carried by this membrane. These inclusions experience the same stochastic fluctuations as 
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those experienced by the membrane itself, resulting in the transformation of their determinis- 
tic (drift) space-time dynamics into a stochastic Langevin-type dynamics described by the Ito 
stochastic calculus. A set of dynamic simulations, resembling the standard MD simulations, 
were performed to investigate the phase behaviour of these inclusions. In addition to stochastic 
forces, these inclusions also experience deterministic interactions described by inter-inclusion 
potentials varying as 1/ R'^ with their separations. The simulation results clearly indicate that 
the capture and aggregation rates change with the temperature and that the embedded mo- 
bile inclusions capture a greater number of the floating inclusions. A further extension of the 
present work would be to include the influence of the surface tension, as well as the bending 
rigidity, by keeping the corresponding term in (j^). This will constrain the fluctuations in the 
surface area of the membrane and would have a direct bearing on the inclusion dynamics. 

The second author(HRS) is grateful to the UK Royal Society for financial support through 
a visiting research fellowship and to the School of Computing and Mathematical Sciences 
(Greenwich University) for their hospitality. Both authors acknowledge useful discussions 
with Professor E. Mansfield FRS on the dynamics of objects supported by surface tension. 
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Figure captions 



Figure 1: Two rod-like embedded (M-type) inclusions vertically placed in an amphiphilic 
fluid membrane. The inclusions impose exponentially decaying thickness-matching constraints 
on the bilayer at the inclusion boundary. Heavy solid lines represent amphiphilic molecules. 
Figure based on 

Figure 2: Two rod-like surface (S-type) inclusions lying on the surface of the membrane. The 
rods have lengths Li and L2, widths ei and €2 and making angles di and 62 with the line joining 



their centres of mass. Figure based on [IT3 . 



Figure 3: A small patch of the membrane showing the stochastic X-Y trajectories obtained 
from equation (|2^) for a sample of four S-type inclusions without the presence of the M-type 
inclusions. Both the drift and diffusion motions can be clearly distinguished. 

Figure 4: A set of snapshots, obtained from dynamic simulations, showing the capture of 
rod-like S-type inclusions (white spheres) at the rod-like M-type inclusion sites (black spheres) 
for static (a to c) and dynamic (d to f) M-type inclusions. The aggregation of the M-type 
inclusions can also be observed (d to f). Only the centres of mass of the inclusions are shown. 

Figure 5: A set of snapshots, obtained from dynamic simulations, showing the capture of 
rod-like S-type inclusions at the rod-like M-type inclusion sites for mobile M-type inclusions 
at T = 100°K (a to c) and T = 350°K (d to f). Only the centres of mass of the inclusions are 
shown. 
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